package scu.maqiang.possion;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

public class PoissionCompositeAsHole {
    public static void main(String[] args) {
        Mesh2T3 mesh2D = new Mesh2T3().square2D(100, 100, T3Type.Left);
        ScalarFunc domainFunc = (xy, label, param) -> {
          double x = xy[0], y = xy[1];
          if (Math.abs(x - 0.5) < 0.25 && Math.abs(y - 0.5) < 0.25) {
              return 1.0;
          }
          return -1.0;
        };
        mesh2D.addDomainAsBoundaries(domainFunc, null, 10);
        mesh2D.toGmsh("Square.msh");
        FES2T31 fes = new FES2T31(mesh2D);
        SRMatrix A = new SRMatrix(fes.GetNdof());
        double[] RHS = new double[fes.GetNdof()];
        double[] x = RHS.clone();
        fes.assembleStiff(new double[]{1.0}, BVPType.COMMON, A);
        fes.assembleSource(new double[]{1.0}, BVPType.COMMON, RHS);
        fes.applyBC_MBN(A, RHS, Direct.All, 0.0, 1, 2, 3, 4, 10);
        NewIterSSolver solver  = new NewIterSSolver(A);
        solver.PCGSSOR(RHS, x, 1.5, 1);
        mesh2D.saveMesh(MeshFormat.msh, "MyPoissonMesh.msh");
        mesh2D.toTecplot("PoissonTwo.dat", x);
        mesh2D.toGmsh("PoissonSolution.msh", x);
    }
}
